Spatially heterogenous dynamics in dense, driven granular flows 
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Interest in the dynamical arrest leading to a fluid solid transition in thermal and athermal 
systems has led to questions about the nature of these transitions. These jamming transitions may 
be dependent on the influence of extended structures on the dynamics of the system. Here we show 
results from a simple driven, dissipative, non-equilibrium system which exhibits dynamical hetero- 
geneities similar to those observed in a supercooled liquid which is a system in thermal equilibrium. 
Observations of the time TR{r) required for a particular particle to move a distance r reveal the exis- 
tence of large-scale correlated dynamical regions with characteristic timescales chosen from a broad 
distribution. The mean squared displacement of ensembles of particles with varying characteristic 
Ti?(r) reveals an intriguing spatially heterogenous mobility. This suggests that a unified framework 
for jamming will have to be based on the connection between the nature of these heterogeneities 
and the effective dynamics. 

PACS numbers: 81.05.Rm, 45.70.-n, 83.10.Pp 



INTRODUCTION 



Materials as diverse as molecular liquids, foams and 
granular matter experience a transition from a fluid-like 
state to a solid-like state characterized only by a sud- 
den arrest of their dynamics. It has been proposed that 
the phenomena associated with this dynamical arrest can 
be explained within the unified framework of a jamming 
phase diagram [T]. Unlike thermal phase transitions, 
jamming does not seem to be characterized by singu- 
larities in any purely static, structural property [21, [J]. 
Numerous experiments and simulations indicate the ap- 
pearance of large-scale dynamical heterogeneities and an 
associated growing dynamical length scale. In colloids 
near the glass transition, for example, fast-moving parti- 
cles were observed to be spatially correlated with a char- 
acteristic cluster size that increased as the glass transi- 
tion was approached [4]. In supercooled liquids, spatial 
inhomogeneities can be identified via a time-dependent 
four-point (two-time, two-space) correlation function [5]. 
Recent experiments on molecular liquids and colloids 
have identified a growing dynamical length scale near 
the glass transition via new multipoint dynamic suscep- 
tibilities 0. In granular materials, spatial structures 
have been directly observed in experiments on flowing 
systems [7, 8]. More recently, measurements of velocity 
correlations in the surface layer of particles flowing down 
an inclined plane have yielded a length scale which ap- 
pears to grow as the flow is arrested @ . While detection 
of such dynamical heterogeneities in numerical studies 
has been more difficult, extended chains of particles ex- 
periencing high frequency collisions are visible in event- 
driven simulations of gravity driven granular flow [10]. 
These structures experience a slow relaxation from colli- 
sional stress at intermediate time scales in a manner anal- 
ogous to temporal relaxations observed in glassy systems 
at low temperatures [11]. There is also some evidence 
for the existence of large-scale structure in images of the 
contact force network in simulations of chute flow [l2ll. 



Several questions naturally arise: (1) Are dynamical het- 
erogeneities a necessary precursor to dynamical arrest of 
the jamming kind? (2) Do these heterogeneities and their 
mesoscopic scale lead to any kind of universal dynamical 
behavior irrespective of significant differences in micro- 
scopic dynamics? (3) Do static structures such as force 
chains observed in jammed granular packings emerge out 
of dynamical heterogeneities? 



In the present work we address these questions via a 
simplified model which focuses on the essential effects 
of driving and dissipation in a flowing granular system. 
Simulations have been performed of a two-dimensional 
gravity driven system of frictionless bidisperse hard disks 
in a hopper geometry. The disks undergo instantaneous, 
inelastic binary collisions and propagate under gravity in 
between collisions; a driven, dissipative system that is 
very different from supercooled molecular liquids. Previ- 
ous work p^] has shown that this minimal model still 
reproduces observable results from related experimen- 
tal systems [l3[. To investigate the relevant time and 
length scales governing the dynamics of the granular sys- 
tem, an interesting analysis for quantifying dynamical 
heterogeneity in non-equilibrium systems [ij] indepen- 
dent of the microscopic interactions between particles is 
employed. Specifically, measurements of the time rR{r) 
required for a particular particle to move a specified dis- 
tance r, provide an appealing physical picture of large- 
scale correlated regions with characteristic timescales 
chosen from a broad distribution. Through a coarsegrain- 
ing scheme first proposed in Ref. ^] which exploits the 
fact that the visible heterogeneity is maximized at a par- 
ticular distance Tc, both the distance Tc (which can be 
associated with average cage size) and the spatial ex- 
tent of the correlated regions ^* are extracted in a com- 
pletely threshold-independent way. Tc can then be used 
to explore the effect of complex collective behaviour on 
dynamic properties such as the average mean squared 
displacement ((Ar)^). 
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II. METHODS 



The grain dynamics used in the simulations are sim- 
ilar to Ref. [15]. Specifically, (i) at each interparticle 
collision, momentum is conserved but the collisions are 
inelastic so the relative normal velocity is reduced by the 
coefficient of restitution (ii) to allow the side walls to 
absorb some vertical momentum we impose the condition 
that collisions with the walls are inelastic with a coeffi- 
cient of restitution jiiwaii in the direction tangential to the 
wall and (iii) since we wish to observe the system over 
many events, particles exiting the system at the bottom 
must be replaced at the top to create uniform, sustained 
flow. Note that collisions are instantaneous and in be- 
tween collisions the particles are driven by gravity. To 
avoid the phenomena of inelastic collapse the coefficients 
of restitution /i and jiiwaii are velocity dependent; if the 
relative normal velocity between particles or between a 
given particle and the wall is less than some cutoff Vcut 
then the collision is presumed to be elastic [16]. The flow 
velocity Vf is controlled by adjusting the width of the 
hopper opening. We also introduce a probability of re- 
flection p at the bottom which reduces the time needed to 
reach steady-state flow. Typically, our simulations were 
done on bidisperse systems (diameter ratio 1:1.2) of 1000 
disks, with ja = 0.8, jiiwaii = 0-5, Vcut = 1 x 10~^ and p 
= 0.5. The simulation was run for a total time of 1 x 10^ 
in simulation time units (smaller particle diameter (ig, 
smaller particle mass and gravitational constant g 
are all set to 1) with the initial time interval of 5 x 10^ 
discarded before recording data to ensure the system has 
reached steady state. During the total time interval of 
500 over which we are evaluating the data, a given par- 
ticle will pass through the hopper 5-10 times depending 
on the flow velocity. 



III. RELAXATION TIMES AND 
CAGE-BREAKING FLUCTUATIONS 

Subdiffusive behavior and caging have become hall- 
marks of the jamming transition, but defining a quan- 
titative measure for these processes which is sensitive to 
the presence of dynamical heterogeneity has been diffi- 
cult. The following observation suggests a possible pro- 
cedure for defining a characteristic time for a particular 
particle which may detect the presence of any spatial het- 
erogeneity in the dynamics: if the relaxation time Tj^{r) 
of a given particle i is defined as the time required for 
particle i to travel a distance r, then any spatial corre- 
lation in r]^(r) (given initial positions of the particles) 
will depend very strongly on the chosen value of r. If r is 
small, then it is not possible to distinguish between those 
particles which are simply rattling back and forth within 
some confined distance and those which may have some 
enhanced motility. On the other hand, if r is too large, 
then a given particle i may have been in regions where 
the local average relaxation time has taken on many val- 
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ues, both fast and slow. In both cases, a spatial map of 
the relaxation time reveals no clear positional correlation 
in regions of r]^(r) beyond some small length scale associ- 
ated with random clumping (Fig.[T^ and c, respectively). 

However, at some intermediate value of r, it should 
be possible to distinguish between localized motion and 
extended translational motion prior to a "washing-out" 
of any spatial correlations between particles with simi- 
lar relaxation times (compare Fig. [TJd, where r = 0.45 
to Fig. [1^ and c, and note the clearer distinction be- 
tween regions with different Tj^{r)). Presumably there is 
some optimal value Tc at which the visible heterogeneity 
in rR{r) is maximized. The primary aim of the analy- 
sis proposed in Ref. [13] is to determine this value of Tc 
and thence the distribution of "critical" relaxation times 
P{t^) where = rR{rc). Particle behaviour such as 
mean squared displacement can then be analyzed in the 
context of for each individual grain and whether or 
not the particle is then "fast" or "slow". 

The method for determining the critical relaxation 
length Tc is as follows: Define the coarsegrained relax- 
ation time T|^(r, /) for a given particle i as the mean local 
relaxation time in a box of size / x / centered on particle 
i. As the size of the coarse-graining region is varied from 
/ = 1 (essentially, the single particle relaxation time) to 
the size of the analysis region / = L, the spatial variation 
in r]^(r, /) should decrease (the result of the coarsegrain- 
ing procedure on the spatial map is depicted for one snap- 
shot of the system in Fig. [2] for coarsegraining lengths of 
/ = 2, 5, 10). Thus, the distribution P{rR{r^ I) should be- 
come increasingly narrower, having zero width at / = L 
where all of the T\^(r^L) = j^^f T^^i^^l) = tr. It is 
clear from the picture that the heterogenities persist up 
to larger / for r ~ 0.4. If the second cumulant m2(r, /) of 
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FIG. 2: Effect of coarsegraining on relaxation time 

map for r — {0.1, 0.4, l.OjSpatial map of the relaxation 
time Tij(r,/) for r = {0.1,0.4,1.0} at Vf = 0.45 for (a) no 
coarsegraining (/ = 1), (b) / = 2, (c) / = 5 and (d) / = 10. 
Colour scale is as for Fig. [1] 



^('?"i?(r, /)) is calculated as: 



m2(r, /) 



{{r},{r.l)-m?) 
((r],(r,l)-r^)2) 



(1) 



then m2(r, /) should be a monotonically decreasing func- 
tion of where the length scale associated with the decay 
of 7712 (r, should vary with r and be maximal at Vc 

It is worth noting at this point that the original mea- 
surement performed in Ref. [3] was for a 2D periodic 
system of soft particles in thermal equilibrium interacting 
via the repulsive potential V{r) = e(cr/r)~^^. While the 
validity of this analysis should not depend on the type of 
microscopic interaction between particles, two additional 
considerations had to be taken into account before per- 
forming this calculation on the flowing granular system. 
Firstly, motion due to the net flow of the system had 
to be subtracted off prior to calculating r]^(r, While 
more complicated schemes involving a spatially varying 
flow field could be devised to more accurately compen- 
sate for this effect, for the analysis here the velocity over 
the region of interest Vf{t) was assumed to be constant 
in space at the snapshot at time t {i.e. so a particle's 
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FIG. 3: Relaxation time distributions P{rR{r,l)) for 
Vf = 1.1 3. Relaxation time distributions P(TR(r,l)) for (a) 
r = 0.4 Tc for this flow velocity) and (b) / = 1 for 
r = {0.1,0.4,0.7,1.0}. In all cases the distributions are for 
Vf = 1.13 and the average relaxation time Tr has been sub- 
tracted off. Note that the average collision time for this flow 



velocity is r ' 

> T. 



1 X 10 and thus typical relaxation times are 



displacement Ar(t) in a given time interval At is calcu- 
lated as Ar(t) = T{t) - T{t - At) - Vf{t - At)At). This 
first approximation seemed to be sufficient to extract Tc- 
Secondly, because relaxation times were only calculated 
for particles in the constant flow region, not all parti- 
cles in the system will be in this region at a given time, 
nor will a given particle necessarily remain in the region 
of interest until it reaches the specified value of r. So 
the calculation only considers a given particle trajectory 
until it departs the analysis region. 

Relaxation time distributions P{rR{rJ)) for r = 0.4 
and / = {1, 2, 5, 10} are shown in Fig.[3K, with a compar- 
ison of the single-particle {i.e., I = 1) relaxation time for 
several values of r in Fig. [3)3. The effect of increasing the 
coarsegraining region for a given value of r is to smooth 
out the distribution for TR{r,l) > Tr, and the width of 
the distribution decreases as expected for increasing /. 
Note that the single particle distribution for / = 1 can be 
quite broad for large values of r. This is a reasonable re- 
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FIG. 4: ^(r) for varying Vf. Decay length scale ^(r) = 
fi 1712(1, r)dl as a function of r for varying Vf. Inset shows 
peak location Tc vs. t'/. 



duced an increase in ^* while Tc remained fixed. Thus 
the crucial length scale which is growing as the flow ve- 
locity decreases in our system appears to be Tc- There 
is some indication of a crossover between a low-velocity 
regime where ^ may be changing more rapidly than Vc to 
a high-velocity regime where changes in Tc dominate as 
described above, but the current data is not sufficient to 
quantitatively distinguish between these two possibilities. 
Additionally, the effect of the finite size of the system on 
the dynamically correlated regions remains an interesting 
unresolved issue. While the correlation lengths measured 
over the range of fiow velocities measured here are still 
sufficiently smaller than the system dimension, certainly 
at the slowest fiow velocity observed the size of these re- 
gions begins to approach the system size (at least in the 
direction perpendicular to the fiow). In order to fully 
observe the behaviour of this length scale as the fiow 
velocity is decreased further, it will be necessary to use 
larger systems. 



suit in the context of the calculation of Tc] the individual 
distributions for a particular value of r 7^ Tc can be broad 
initially, but unless there is any significant spatial corre- 
lation in relaxation times the spatial structure of ri^(r, /) 
will vanish rapidly with increasing / compared to a value 
of r with nontrivial spatial correlations. Thus the utility 
of TR{r^l) lies in its ability to emphasize spatial corre- 
lations between fast and slow particles {i.e. if there are 
particles of widely different mobilities in the system but 
they were randomly distributed in space then no inter- 
esting correlation length is visible). 

The length scale associated with the decay of 7x12(1 ^r) 

can be calculated from ^(r) = m2(l^r)dl. As can be 
clearly seen in Fig. [H there is a distinct peak in ^(r) 
from which two separate and interesting length scales 
can be extracted . One is the location in r of the peak; 
this quantity is Vc and it is the length scale over which 
the dynamical heterogeneity in relaxation times is max- 
imized as described above. An additional interpretation 
can be made of Vc as the critical fiuctuation associated 
with cage-breaking in supercooled fiuids (or equivalently, 
the typical cage size). This was proposed in Ref. [14] and 
confirmed experimentally for slow, quasistatic granular 
shear fiows in Ref. [l7]. Referring to the inset to Fig. HJ 
this quantity is decreasing monotonically with fiow ve- 
locity, although the exact functional form is difficult to 
determine. It is also important to note that because Vc 
can be related back to a critical relaxation time, changes 
in Tc with Vf will correspond to changes in the distribu- 
tion of r^. 

The second length scale which may be determined from 
^(r) is the peak height which represents the average 
size of the heterogeneity From Fig. [H it is difficult to 
determine if ^* is varying with Vf^oi varying very slowly 
compared to the changes in Vc This is in contrast to 
the results of Ref. [14,] where an increase in density pro- 



IV. SPATIALLY HETEROGENOUS MOBILITY 

With the critical fiuctuation Vc determined, the crit- 
ical relaxation time for each particle can be defined 
(Fig. [1)3 shows a snapshot of the Vf = 0.45 system at 
r = Tc). As in Ref. [14], different particles can now be 
grouped into subsets according to their value of and 
the contribution of these subsets to dynamical measures 
(in this case, the mean squared displacement (MSD) 
(Ar^)) can be evaluated. Plots of (Ar^) vs. time are 
shown in Fig. [5] for varying fiow rates (the distribu- 
tions of critical relaxation times for a given fiow rate are 
shown in the insets to Fig. [5]). In all cases, the MSD 
averaged over all particles in the system in the short- 
time regime is approximately ballistic ((Ar^) ~ t^'^). 
The long time power-law behaviour is characterized by a 
fiow velocity dependent exponent ((Ar^) ~ t^('^f\ where 
0.9 < ^(vf) < 1.2 in the range of fiow velocities studied) 
which is nearly diffusive. Fig. [5^ shows the power-law 
fits to the mean squared displacement curves for two dif- 
ferent fiow velocities. At times comparable to the aver- 
age time it takes a particle to cross the analysis region 
tcross = l/vf^ the particle falls out of the dense constant 
fiow velocity region and accelerates under gravity, leading 
to the sharp increase in the mean squared displacement 

at t ~ ^cross' 

The MSD curves for "slow" (r^ ~ t cross) and "fast" 
(r^ ~ (1/3 — 1/2) t cross) particles reveal a deeper com- 
plexity not visible in the mean squared displacement av- 
eraged over the system as a whole. As '^j ^ 0, these two 
groups of particles show increasing separation from the 
average MSD, with the fast particles experiencing greater 
displacement than the average in a given time interval, 
and the slow particles exhibiting a plateau (Fig. [5)3-c). 
For the slowest fiow velocity measured, where the plateau 
is best defined, the plateau occurs at (Ar^) ^ r1 , corre- 
sponding to the experimental result of Ref. [13 where 
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FIG. 5: Mean squared displacement (Ar^) for varying 
Vf, (a) Mean squared displacement (MSD) (Ar^) (() indi- 
cates average over all particles in a snapshot at time to, then 
average over all snapshots to) for Vf — 2.05 and Vf — 0.45. 
(b)-(c) Comparison of "slow" particles (r^ ^ t cross) and 
"fast" particles (r^ ^ 1/3 — 1.2tcross) to the average over all 
particles (c/. discussion in text) for (b) Vf — 1.13 and (c) 
Vf = 0.45. 



■ ha8 |been directly associated with the cage size. Thus 
transport in this system is governed by the 
heterogeneity: fast and slow regions are ev- 
depicted pictorially in Fig. [T] with particles in 
region experiencing caging. As the flow velocity 
, and the length scale Vc increases, the effect of 
^esoscopic regions becomes more pronounced, with 
w particles at = 0.45 remaining caged for the 
of their time to flow through the system (recall 
all of these displacement measurements motion 
the average flow has been subtracted off). This 
of rich collective dynamics occuring even in the 
abserl(0© of any evident density fluctuations was also sug- 
gested by Menon and Durian in experiments on dense 
:-driven granular flow using diffusing wave spec- 
)y (DWS) [18]. This is also in analogy to the MSD 
Lennard- Jones system of Ref. [14], where slow 
s (defined as the particles in the slowest 40% of 
were seen to be caged over some time Tmix dur- 
l3ich the fast particles experienced diffusion with a 
i(>n constant twice that of the assembly of disks as 



CONCLUSIONS 



t (Simulation Time Units)^ summary, an interesting analysis for quantifying dy- 
namical heterogeneity in non-equilibrium systems [14] in- 
dependent of the microscopic interactions between par- 
ticles Jias proven to be a useful tool for further inves- 
into the relevant time and length scales gov- 
the dynamics of the gravity-driven granular sys- 
CVIeasurements of relaxation time ri^(r), defined 
time required for a particular particle to move a 
[ed distance r, provide an appealing physical pic- 
f large-scale correlated regions with characteris- 
escales chosen from a broad distribution of re- 
times P{tr). Through a clever coarsegraining 
which exploits the fact that the visible hetero- 
imnffiit}/! is maximized at a particular distance Tc, both 
^?he di^^ance Vc (which can be associated with average 
cage size) and the spatial extent of the correlated regions 
^* were extracted in a completely threshold-independent 
way. While ^* was shown to be ~ 5 — Od independent 
of flow velocity, the cage size Vc grows as the flow ve- 
locity decreases. Vc then specifically defines the critical 
distribution of relaxation times P{rf^) which can be used 
to explore the effect of complex collective behaviour on 
dynamic properties such as the average mean squared 
displacement ((Ar)^). The MSD in turn reveals a non- 
trivial spatial dependence in the mobility, indicating that 
while the total mean squared displacement at long times 
may seem approximately diffusive, there is underlying 
complexity in the form of "fast" and "slow" regions in 
the sample. This result is in agreement with several 
experimental studies on dense, driven granular materi- 
als [l7|, ll8[ and other jamming systems [4]. Addition- 
ally, there is a qualitative resemblance between the spa- 
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tial maps of relaxation time near the critical value 
{c.f. Fig. [1)3) and recent measures of dynamic propen- 
sity and local Debye- Waller factors in 2D simulations of 
glass- formers [19]. 

One as yet unresolved puzzle is the relationship between 
the frequent-collision/high-stress chains identified in pre- 
vious studies and the mesoscopic clumps of fast 
and slow particles characterized by a particular value of 
r^. From observations of the spatial map of collision fre- 
quency and critical relaxation time at a given snapshot of 
the system (not shown) it is difficult to ascertain if there 
is any correlation between the chains and regions of vary- 
ing r^, indeed, there can still be strong heterogeneity in 
the critical relaxation time in the absence of any chains. 
One measureable connection is that (r^) ~ 1/'^/, which 
is the same scaling with ffow velocity that was measured 



in the stress autocorrelations in Ref. [ll| and associated 
with the lifetime of the chains. Thus TRVf = constant 
could be an important clue to the dynamical principles 
leading to the heterogeneities, and may provide some ad- 
ditional insight into the reasons for why the frequently- 
colliding/high-stress chains form in these dense granular 
systems. 
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